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Abstract 

We extend the generating function technique for calculation of single molecule photon emission 
statistics [Y. Zheng and F. L. H. Brown, Phys. Rev. Lett., 90,238305 (2003)] to systems governed by 
multi-level quantum dynamics. This opens up the possibility to study phenomena that are outside 
the realm of purely stochastic and mixed quantum-stochastic models. In particular, the present 
methodology allows for calculation of photon statistics that are spectrally resolved and subject 
to quantum coherence. Several model calculations illustrate the generality of the technique and 
highlight quantitative and qualitative differences between quantum mechanical models and related 
stochastic approximations. Calculations suggest that studying photon statistics as a function of 
photon frequency has the potential to reveal more about system dynamics than the usual broadband 
detection schemes. 
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I. INTRODUCTION 



Single molecule spectroscopy (SMS) has become a versatile and powerful tool for the 
study of condensed phase systems in chemistry physics and biology [1-12]. Unfortunately 
the very qualities that make SMS such a powerful technique, have also led to significant 
theoretical challenges in describing experimental data. The ultra-microscopic nature of the 
physical systems under study leads to randomness in the behavior of experimental signals 
due both to thermal agitation of the photoactive portion of the system and the inherent 
randomness of the spontaneous emission process itself. While SMS has been hailed for its 
ability to probe these fluctuations directly, it remains difficult to extract physical pictures 
for molecular dynamics based solely on SMS data streams. Some of this difficulty is likely 
fundamental (current SMS experiments may not collect sufficient data to allow for direct 
inversion to molecular dynamics), but even if SMS data were sufficient to differentiate be- 
tween all viable physical hypotheses, it remains an open question as to the best means to 
simulate such models to allow for comparison with experiment. Indeed, much effort has 
been expended on the theory of interpreting/modeling SMS trajectories, particularly in the 
context of stochastic models for chromophore dynamics [12-30]. Stochastic models, though 
certainly illustrative and powerful, ultimately face certain limitations in the modeling of 
phenomena that are inherently quantum mechanical, such as spectroscopy. Quantum co- 
herence can not be captured, quantization of nuclear eigenstates is not naturally formulated 
within a stochastic scheme, and the parameters of stochastic models are often difficult to 
equate with their microscopic origins. As the following work will show, even stochastic mod- 
els systematically derived from underlying quantum considerations can lead to quantitative 
and qualitative differences from fully quantum calculations. 

Until recently, Monte Carlo Wave Function simulations (MCWF) [31] and related tech- 
niques provided the only fully satisfactory route to theoretical calculation of single molecule 
photon counting observables [32] including quantum mechanical effects. A few other studies 
have touched on specific aspects of quantum dynamics applicable to SMS [24, 27, 33-35], 
but without complete generality. Recent work by us [36-39] and others [40-44], has estab- 
lished generating function techniques as a general means for calculating statistical quantities 
of single molecule photon counting experiments. The only fundamental limitations to this 
approach are that you must consider the spontaneous emission of photons to be governed by 
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rate processes and the directly calculated quantities are statistical moments of the number 
of photons emitted [36, 37, 41]. 

The bulk of previous work with the generating function approach has focused on two level 
chromophores with stochastic modulation by the environment, however the method is equally 
applicable to multi-state quantum systems. The extension to multi-state quantum systems 
was suggested by us [37] and formally carried out by Mukamel [41]. Sanda and Mukamel 
[44] have recently used the generating function approach to derive formal perturbative ex- 
pressions (in the applied field strength) for low order photon counting moments. Though 
interesting from a theoretical standpoint, the derived expressions are complex enough that 
implementation will be impossible for all but the simplest model systems (second order mo- 
ments require solution of a six point quantum correlation function, higher moments need 
larger correlations). As a numerical technique, the generating function approach has promise 
to study varied systems without limitation to low field strengths. 

The present paper considers several model systems to demonstrate the use of the generat- 
ing function approach as a numerical tool for predicting SMS photon counting observables. 
In addition to calculation of photon counting moments for broadband detection schemes, 
as has been considered previously, we also calculate emission statistics for photons specific 
to particular molecular transitions and degenerate sets of transitions. For systems where 
vibrational structure is well resolved compared to natural line widths, this is equivalent to 
the calculation of spectrally resolved emission statistics. From a conceptual and numerical 
standpoint these calculations are no more difficult than broadband detection calculations. 
The simulations we have carried out suggest that significantly more information stands to 
be learned from photon counting experiments when photon statistics are broken down by 
color. 

This paper is organized as follows. Section II presents the underlying theory and notation 
necessary to introduce our calculations. Although there are many details to be considered 
here, the conceptual framework for calculating photon statistics in the many-level case is no 
more complex than for two level chromophores. Given the reduced Liouvillian operator for 
density matrix dynamics of the chromophore system, calculation of the generating function 
for photon number and/or low order statistical moments is immediate. Most of sec. II is 
dedicated to describing the Liouvillian operator itself, not the extension of this matrix to cal- 
culation of experimental observables. Sections III and IV present numerical calculations for 
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chromophores coupled to a two level system and an harmonic vibrational coordinate. Many 
different regimes are considered, both to display the flexibility of the present formulation in 
numerical calculations and to highlight differences between fully quantum calculations and 
commonly employed stochastic approximations. In sec. V we conclude. 

II. THEORETICAL BACKGROUND 

A. General considerations for chromophore dynamics 

The picture we present is the natural extension of the optical Bloch equations to multi- 
level quantum systems in a condensed phase. Our methodology has been adopted both to 
make connection with our previous work on two level chromophores [36-38] and because the 
necessary theoretical/computational tools for chromophore dynamics are well established in 
the literature. 

We imagine a single chromophore in a condensed phase environment driven by an external 
laser field. It is assumed that the field is strong enough to warrant a classical treatment of 
this perturbation so that dynamics, in the absence of any other system-field interactions, 
would be dictated by 

m = - l T [H sy \p\ + l - h E(t)-[^p\. (i) 

In the above H sys is the Hamiltonian for the unperturbed chromophore-environment system, 
fi is the electric dipole moment for this system and E(t) is the classical applied laser field. 
pit) specifies the density matrix for the molecule only. This evolution assumes the radiation 
is of sufficiently long wavelength (and the chromophore sufficiently localized) to allow the 
dipole approximation. 

What the above dynamics neglects is the relaxation of the driven molecular system. The 
coupling between system and the quantum radiation field provides a route for this relaxation 
to occur: the spontaneous emission of photons. It is these photons that are registered in 
SMS experiments and hence inclusion of the spontaneous emission process is absolutely 
essential. Within the standard approximations, the quantum radiation field is integrated 
over to provide rate constants for emission of photons between various molecular transitions 
[45, 46]. This leads to a master equation approach for incorporating emission events as pure 
rate processes. The rate for spontaneous emission of a photon, causing a jump from system 
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eigenstate i to eigenstate j, is calculated by application of Fermi's golden rule (using the 
coupling between system and quantum radiation field as the perturbation) 




(2) 



D a = (*IAU>- 



The collection of constants appearing in this expression have their usual meaning, but we 
will not be concerned with them in this work. What is important to us is the dependence 
on the transition dipole moment Dij, which serves to mediate relative rates of emission for 
different chromophore transitions. In principle, energy splittings (u)ij) impact the rates as 
well, but we shall be concerned with electronic transitions where differences in this quantity 
between various allowed transitions are much smaller than the splitting itself. In this limit 
we expect inconsequential variations on the basis of energy differences. 

Perturbation theory applied to the entire system density matrix evolution (as opposed 
to just a single rate calculation) additionally tells us that the population lost from state 
i, via Tij decay, ends up in state j. Also, it specifies that the i — > j transition causes all 
associated coherences (pik, pki) to decay at the rate of 1^/2. The net effect of all spontaneous 
emission processes in the system is the additive contribution of these three effects (loss of 
population from state i, gain in population of state j and loss of coherence for all allowed 
% — > j transitions.) We neglect radiative level shifts in the system states and ignore all other 
couplings (virtual photon transitions) caused by the presence of the quantum radiation field. 
These other couplings are unimportant when system energy levels are non-degenerate as the 
implied perturbations are non-secular [45, 46]. The non- degeneracy condition is met by the 
systems studied in this work. 

Keeping those contributions specified in the last two paragraphs, implies that we sup- 
plement our chromophore equations of motion with non-Hamiltonian evolution terms corre- 
sponding to spontaneous emission. The form of this augmentation is most transparent in the 
basis of H sys eigenstates. Rewriting eq. 1 in this form yields (The summation convention 
over repeated density matrix labels is assumed throughout this work.) 



mPm + L i:j . kl (t)pki + h^pki. 



(3) 
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Here, h sys and L, E (t) are Liouville super-operators (matrices) corresponding to the commu- 
tator expressions in eq. 1 [47]. Note that our definition incorporates the factor —i/h within 
h sys and h E (t). L r is the matrix effecting spontaneous emission processes. The elements of 
this matrix are provided by the arguments of the preceding paragraphs (i ^ j assumed in 
the following) 



with all other elements zero. 

In what follows, it will be convenient to partition the matrix L r into its positive and 
negative pieces, so that 



with L +r consisting of the terms specified by the second line of eq. 4 and L _r comprised of 
the remaining terms from the first and third lines. 

One final important point is that while eq. 3 provides effective dynamics for the system 
with implications of field fluctuations handled implicitly, this dynamics will still be far too 
complicated for exact practical treatment when the system of interest is composed of a 
chromophore embedded in a condensed phase. The problem is simply one of a complex 
dynamics associated with a quantum mechanical many body system. When it is possible 
to make some effective separation between the relevant part of the system and a weakly 
coupled (and fast) bath this problem can be overcome in exactly the same method employed 
to remove the radiative field from explicit consideration. Writing 



for a "system" Hamiltonian composed of two parts: ch (the chromophore which is directly 
coupled to the applied field) and b (the bath) weakly coupled by V, we arrive at an equation 
of motion for the reduced chromophore density matrix a through application of standard 




(4) 




L r = L +r + L~ r 



(5) 



H sys = H ch + H b + V 



(6) 
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Redfield theory [48, 49] 



&ij(t) = -i^ijCTij + L^. H (t)(TM + L^-. w CTfe/ + Rij-klCTkl = ^ij;ki(t)cr k i. (7) 

Here, the matrix R is the usual Redfield matrix to account for bath perturbations on the 
chromophore and the matrix h(t) reflects the entire dynamics for a. We note that additivity 
of contributions stemming from quantum field, bath and classical (laser) field perturbations 
to the dynamics of the chromophore should be viewed as an approximation of "independent 
rates of variation" [46]. We neglect frequency shifts of the chromophore due to V, so that 
the labels ij now correspond to eigenstates of H ch . We consider this set of approximations 
as the natural extension of the optical Bloch equations to multi-level systems in a condensed 
phase. Specification of the matrices h E , L r and R will allow us to apply this formalism 
to various physical problems and several model systems will be considered in the following 
sections. 

B. Extraction of photon counting moments 

Extending the picture of the preceding section to calculation of photon counting statistics 
for single molecule measurements proceeds in a manner analogous to the case for two level 
chromophores [37, 40]. The formal solution has been presented in ref. [41] and we present 
here a brief derivation following ref. [37] to clarify our notation and to extend this picture to 
the calculation of photon counting moments for individual spontaneous emission transitions 
(as will be useful in spectrally resolved emission spectroscopy). 

Imagine a detector capable of differentiating between photons that are emitted for par- 
ticular chromophore transitions. In certain cases this would be possible by only selecting 
photons within a certain frequency window, in other cases this might not be experimentally 
feasible but should be regarded as a gedanken experiment. That portion of L(t) responsible 
for placing the chromophore in a lower energy state immediately following the transition 
of interest is of special importance for calculating statistics associated with this transition. 
From eq. 4 this is the element L^ r aa with the numerical value T a b, assuming that we are 
following a — > b emissions. Partition eq. 7 to give this single part of the evolution a unique 
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status 

&ij(t) = L' ij . kl (t)a k i + T ab 5i j)bb 5 k i- aa <7 M = \J i:j . kl {t)cr kl + L+.^f o kl (8) 

where L'(t) is that portion of L(t) not pulled out in h +Fab . In exact analogy to the case 
with only a two level chromophore, it is the operator L +ra6 that dictates when an a — > b 
spontaneous emission event occurs. Following exactly the same arguments as in ref. [37] 
allows us to write 

where is that portion of the reduced density matrix corresponding to systems that have 
previously emitted exactly n photons via a — > b transitions. 

To facilitate the extraction of photon counting moments, we introduce a generating func- 
tion version of eq. 9 

g^t, s) = u trM (t)g kl (t, s ) + sL+Ftfg^t, s) = Li^t, s)g kl (t, s ) (10) 

oo 

g(t,s) = J2^ [n) (t). 

n=0 

The actual generating function for a — > b photon emissions is obtained by summing over all 
"population" elements of g(s,t) 

G(s,t) = Y t g ii (s,t) (ii) 

i 

which allows for the usual extraction of probabilities for n photon emissions [50] 



1 d n 



(12) 



s=0 



and factorial moments [50] 



am 

(nW)(t) = (n(n - \){n - 2) . . . (n - m + l))(t) = ^G(s, t) 



(13) 



s=l 



Our primary concern in this work shall be the calculation of moments. To this end, we 
differentiate eq. 10 with respect to s yielding equations for the d m /ds m g elements which, 
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when summed over population elements, yield the moments (when s = 1). 



d fd m g(s,t)\ w ,fd m g\ V . T d m - x g 

-7T — — - = h(t, s) — — + mL +r ° b - ^ 

dt V ds m J V ds m / V ds" 1 " 1 



(14) 



The high order derivatives are dependent upon all lower derivatives as can be seen by iter- 
ating this equation. For example, moments up to and including second order are generated 
by solving the set of equations 



dQ{s,t) 

ds 
d 2 G(s,t) 



( 



L(i, s) 
L+ r " 6 (s) h(t,s) 
2L +r «»(s) h(t,s) 



ag(s,t) 

ds 

a 2 G(s,t) 



(15) 



Evaluation at s = 1 provides the moments up to second order by way of eq. 13. Since L(t, s) 
and L +ra6 are N 2 x iV 2 matrices for a quantum system with N states, the above expression 
corresponds to solving 3iV 2 coupled equations. In the cases considered in this work, we will 
take E(t) to have sinusoidal time dependence so that the explicit time dependence within 
L(i) may be removed by moving to a rotating reference frame and applying the rotating 
wave approximation (RWA). In this case, solution of these equations is easily accomplished 
by directly exponentiating the 3iV 2 x 3iV 2 matrix as outlined in the next section. Equation 
15 is central to all results in this paper and, in principle, could have been directly solved to 
reproduce all the calculations presented below. In practice, we used a numerically simpler 
scheme to obtain our results derived from eq. 15. This numerical technique is elaborated on 
in sec. II D. Formation of the matrices L(i, s) and h +Vab for use in any numerical scheme 
follow from the preceding section. Specific choices for these matrices depend upon the 
physical systems under consideration and will be detailed with presentation of our chosen 
applications. 

The above derivation has assumed that we are interested in the statistics of photons 
emitted from one particular chromophore transition (a — > b). When we are interested in 
broadband detection with all photons counted equivalently, the structure of eq. 15 remains 
unchanged. However, the matrices L(t, s) and L +rab have different forms. In that case we 
substitute L +r for L +ra6 and L(t, s) is now the matrix formed by appending s to every spon- 
taneous emission matrix element within L(t) having a positive sign (i.e. the whole of L +r ). 
Calculation of moments for photons associated with some subset of transitions (perhaps 



transitions inside a certain frequency window) proceeds by generalizing to placement of s 
variables only on the elements associated with the relevant transitions and making the cor- 
responding changes to L +r . In principle, we could introduce a number of different auxiliary 
variables - each variable corresponding to a particular transition or subset of transitions. 
This leads to expressions for cross correlations between various transitions. The extension 
is straightforward, but not explicitly presented here as we do not calculate any such cross 
correlations in this work. 



C. Model Hamiltonians and practical considerations 

In this work we shall be concerned exclusively with model systems consisting of a chro- 
mophore with two electronic states (ground \g) and excited |e)), so that 



H ch =\g)H g (g\ + \e)H e (e\- (16) 



H g and H e are, respectively, the chromophore Hamiltonians for nuclear motion within the 
ground and excited states, with eigenfunctions and eigenvalues specified by 

H g \n g ) = e ng \n g ), (17) 
H e \m e ) = e m Jm e ), 



for m e — 1 . . . N e , rig — 1, . . . N g . In our numerical applications, we consider only a finite 
number of eigenstates associated with nuclear motion, and adopt the convention here. The 
nuclear ground state in the excited manifold is assumed to lie higher in energy than the 
nuclear ground state of the ground manifold by an amount hu eg . It is to be understood 
that this chromophore Hamiltonian dictates dynamics in the sense implied by eq. 6. H ch 
is responsible for the evolution that we designate to be the most important to chromophore 
dynamics. The effect of the environment (bath) will be felt through coupling dictated by V. 

Interactions with the radiation field depend upon the matrix elements of the system's 
dipole moment operator as evidenced by eq. 2 and the presence of ft in h E (t). We treat 
these matrix elements in the Condon approximation [47] such that 



D 



n g ;m e 



(g\fi\e)(n g \m e ) = /j, {n g \m e ). 



(18) 
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The dipole operator is assumed to act solely in the electronic space with only off-diagonal 
coupling between ground and excited states. Individual transition intensities are mediated 
by the overlap of nuclear wavefunctions. We always consider a monochromatic exciting field 
of constant intensity and polarization direction, so that 

E(t) = S cos(iu L t). (19) 

For future notational simplicity we define constants T and flo as 

(20) 



37re h(? 
Q = £ - /j, /H 

These constants represent the spontaneous emission rate and Rabi frequency for an electronic 
transition between states with perfect overlap of nuclear wavefunctions. 

These definitions allow us specify the form of matrices L E (t) and L r . L r follows imme- 
diately from eq. 4. All we need are the emission rates for all % — > j transitions. Since 
our models only allow transitions between excited and ground electronic states we need only 
consider rates of the form r| e >| me ) ; | 3 }|„ 9 ) = T merig with values 

r me „ 9 = T \{m e \n g )\ 2 . (21) 

All positions in L r diagonal in the electronic subspace are necessarily zero due to our as- 
sumptions about the dipole operator, so the above completely specifies the L r matrix. 

Formation of ~h E (t) is slightly more complicated due to the nature of the coupling to the 
applied field, which makes for a matrix less sparse than the emission matrix. We first realize 
that, as in the usual optical Bloch equations, density matrix elements diagonal in the elec- 
tronic subspace are coupled to those off-diagonal in the electronic subspace and vice versa. 
Also, by analogy to the optical Bloch equations we retain only those terms corresponding 
to resonant excitation by the field (i.e. a photon is absorbed and electronic state rises or 
a photon is emitted and state drops) by invoking the Rotating wave approximation (RWA) 
[46] . We make use of the definition 

tt me n g = £l \(m e \n g )\ (22) 



11 



to give the elements of h E (t) within the RWA 



h E = 


-Lf , ■ 


+ 2^™9 fc e e L 3m g ,lg 


(23) 


h E . = 


fcgle jTlgYTlg 


q —tuj^tir 

2 ' em 9 n g>kg 




h E . = 




= ~2^l g m e C L 3n e ,k e 




h E . = 


-LP 







(24) 



Over bars represent complex conjugation. (The above definitions assume that our dipole op- 
erator matrix elements are real quantities. In the presence of a magnetic field this condition 
could be violated, but we restrict attention away from such cases.) 

The only portion of h(t) remaining to be specified is the Redfield matrix for transitions 
of the chromophore induced by environmental bath fluctuations, R. The relaxation matrix 
is given by [48, 49] 

%;fcZ = Sik ^2 ~ ^3 ^2 ^rrk + t jlik + tfliki ( 25 ) 

r r 

where 

i r°° . „ „ 

*S* = ^ J dre-^ <Vy(r)V tt (0)> 6 , (26) 

are Fourier-Laplace transforms of the correlation functions of the system and bath coupling 
at the specified frequency. The bath-space Heisenberg operators are defined by 

V ki (r) = ei^Vke-i"\ (27) 
Vu = (k\V\i) 

and the averages (. . specify a thermal average over bath degrees of freedom only. In 
all models we consider, bath fluctuations are capable of causing transitions between levels 
within a particular electronic state, but are not permitted to induce radiationless transitions 
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between electronic states. Further discussion on the evaluation of K. will appear in sections 
III and IV as specific models for chromophore and bath are introduced. 

Given the particular form of our model systems, it is highly beneficial to solve eq. 7 in a 
rotating reference frame by introducing new variables 



'n g m e 



'n e m g 



&n g m, e & 



(28) 



(T P 



UJ£,t 



a 



n e m e 



<Jn a m a ■ 



The primary advantage of this formulation being that eq. 7 is recast in a form without 
explicit time dependence 



&ij(t) — —^ij;kl Cr kl + ^ij;kl a kl + ~^ij-kl a kl + 7Hij;kl a kl = Uj;kl&kl 



(29) 



where the diagonal matrix W is given by 



w 

" n g m g ;n g m g 


— Un g m g 




W 

vv n e m e ;n e m e 






W 


^n e in g 


~ U L 


W 

n g m e ;n g m e 


= ^n g m e 


+ OJL 



(30) 



The matrix h E is simply the matrix specified by eq. 23, evaluated at t — and the remaining 
matrices are unchanged relative to the original basis. Since the populations of a are identical 
to a, we may calculate photon emission statistics using these transformed variables without 
any changes to the formalism of the preceding subsection. In particular, we may calculate 
eq. 15 as 



/ 



dd(s,t) 



\ 



ds 
d 2 d(s,t) 
\^s^ ) 



( 



Us) 



L+ r ° 





L(s) 








^ 2L +r ^(s) L(s) J 



\ 1 Q(s,t)^ 

dg(s,t) 

ds 

a 2 g(s,t) 

\ ds 2 ) 



(31) 
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where the time independent L is specified by eq. 29 and G(s,t) is given by 



oo 



(32) 



n=0 



Summing over the "population" elements of Q still returns the original generating function 
for photon emissions, G(s,t), so calculations in this frame return emission statistics equiv- 
alent to the original formulation. Numerics in this basis are preferred, since eq. 31 may be 
solved simply by direct matrix exponentiation. 

D. Reported quantities and numerical details 

The bulk of the preceding sections has been devoted to establishing models for reduced 
chromophore dynamics, i.e. how to specify the superoperator matrix L(t) in eq. 7 or 
the corresponding time-independent matrix L in eq. 29. Given this matrix, it is a trivial 
programming task to extend the standard calculation of density matrix evolution to photon 
counting observables. The matrix L(s) is formed by appending the auxiliary variable s to 
elements of L +r reflecting spontaneous emission transitions of interest. In the case of a 
single relevant transition, only one element is modified. In broadband detection we append 
an s to the entire L +r matrix. Given L(s), the block form of eq. 31 follows immediately 
and calculation of Q is provided by simple matrix exponentiation. Summing over population 
elements of d m Q(s,t) / ds m for s = 1 yields the factorial photon counting moment of order 
m. Although the matrix in 31 is specific to calculation of m — 2, higher order moments 
can be calculated in analogous fashion by extending the block matrix as implied by eq. 15. 
Since we assume no photon emissions prior to t — 0, the initial condition employed in eq. 
31 is simply Qij(s, 0) = £7^(0) with all s derivatives of Q equal to zero. 

The moments reported in this work will be presented in terms of absorption and emission 
lineshapes and corresponding Mandel Q parameter [51] spectra. Mandel's Q parameter is 
related to the factorial moments via 



and serves as a convenient means to report second order photon statistics. Positive Q values 



Q(t) 



(n 2 )(t)-(n) 2 (t) 
(n)(t) 



-1, 



(33) 
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reflect photon bunching behavior (an elevated variance in n relative to Poisson processes 
with the same mean), negative Q values anti-bunching behavior (diminished variance in n 
relative to a Poisson process with the same mean) and Q = is consistent with purely 
Poissonian statistics. 

Energy conservation implies that we may calculate absorption lineshapes, by counting 
the relative rate of photon emission (photons from all transitions are counted) as a function 
of the exciting frequency 



Every emitted photon corresponds to a prior excitation of the chromophore, and hence a 
quantum of energy (hu^) extracted from the incident field. We evaluate lineshapes in the 
limit of long times to insure that the system is in a steady state. The time dependence of 
d{n)/dt at early times is interesting as well [36, 37], but not specifically considered in this 
work. The Q parameter absorption spectra are calculated in analogous fashion, although 
the definition of Q, with (n)(t) in the denominator, insures saturation to a constant value 
as time becomes large. It is unnecessary to take a time derivative to report a meaningful 
quantity here and the Q parameter itself as a function of exciting frequency is reported. 
Again, in the "absorption Q spectra" we collect all photon emissions (broadband detection). 

Emission lineshapes and Q parameter are calculated in similar fashion, but we resolve 
the photon statistics by frequency of the emitted photons. More precisely, we resolve by the 
transitions the photons originate from. In the cases we consider, the allowed transitions are 
either well resolved in frequency (frequency differences much larger than natural linewidths) 
or perfectly degenerate, so that there is no ambiguity in assigning photons to a particular 
frequency "window" . We report our results as 



The above notation specifies that we only consider photons from transitions on resonance 
with uo E . Collection of these statistics follows the prescription previously described. The 
matrix L(s) depends on cue as placement of s variables is dictated by which transitions 
are on resonance with u E . We note that our emission "spectra" are thus not quite spectra 




(34) 




(35) 
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in the usual sense. Our spectral lines are infinitely sharp, without broadening (see fig. 
5). In principle, we could artificially broaden these lines by making them Lorentzians with 
the natural linewidth of each transition, but we have not done so. What our calculations 
directly provide are the statistics associated with particular molecular transitions, not the 
actual frequency of the emitted photons. Note that our lineshapes will also, in general, 
depend upon the frequency of the exciting light as different excitations can lead to different 
steady state populations of the chromophore. 

The Q parameter emission spectrum follows similarly 

r\i \ — ^^i= u «^ — ( n uij=ujE) 2 , s 

Q(uj e ;lu l ) = r 1 (36) 

where we stress that the photon numbers n collected above reflect only those photons stem- 
ming from transitions on resonance with u>e- 

For multi level quantum systems the matrix of eq. 31 can become very large (3iV 2 x 3N 2 
for N quantum levels). If moments higher than second order are desired, the matrix becomes 
even bigger. Direct exponentiation of such matrices over a wide range of frequencies is 
computationally expensive and, for sufficiently large N and/or moment order, eventually 
becomes computationally intractable. In this work we focus on statistics calculated in the 
long time (steady state) limit. For direct exponentiation, this limit has the additional 
computational complications associated with the identification of a time sufficiently large 
for the steady state to be attained, yet sufficiently small to insure numerical stability. When 
only steady state information is desired, analytical progress can be made on eq. 31, allowing 
calculation to proceed via diagonalization of matrices no larger than N 2 x N 2 and without 
the need to identify a suitable finite time at which the long time limit is reached. The 
calculation is summarized below. 

The equations of motion for Q and its s derivatives (Eq. 31) can be formally integrated 
to yield 

Q\t) = f 'dt'e^h +r p s . s . (37) 
Jo 

g"(t) = 2 f dt'e h ^h +r f df'e^'-^h^p^. 
Jo Jo 

Here we have assumed that the system began in the steady state at t — and that we began 
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counting photons at t — (different initial conditions lead to negligible corrections in the 
long time limit). We have introduced a prime notation for s derivatives (i.e. ^ = Q') and 
we have evaluated everything for s = 1. The steady state limit for the density matrix p s . s . 
is expected on physical grounds for systems driven by external perturbations and allowed 
to relax via radiative and non-radiative transitions - its existence was verified for the model 
systems studied in this work. 

The matrix L may be diagonalized and we write A = x~ x Lx with A the diagonal repre- 
sentation of L. The columns of x consist of the right eigenvectors of L and the rows of x~ l 
are the left eigenvectors of L. The associated eigenvalues of L are complex numbers with 
negative real parts, excepting the single eigenvalue associated with the steady state which 
is zero. Ordering the eigenvalues {\ s . s = 0, A2, A3, . . .}, so that 



A 




A 2 
A 3 



\ 



(38) 



we see that it is possible to partition the time evolution operator U(r) = e Lr = U + U\(t) 
into two pieces such that the first corresponds to the (lack of) evolution of the steady state 
and the second piece reflects all other dynamics in the system. 



U = x 



1 





\ 



Ui(t)= X 



^00 ...^ 

e A2T ••• 
e A3T ••• 



V 



X 



J 



Partitioning the matrices in this way allows us to explicitly carry out the integrations in 
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eq. 37 to give (large time limit assumed) 



Q' = (tf/ + X)L +r ( 

G" = t 2 (U h+y Ps . s 

f 



Ps.s. 




(39) 
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The long time (steady state) limit for the rate of photon emission (intensity) and the Q 
parameter follow immediately 



where the summations are over the population elements of the resulting vectors. 

Eq. 40 was used in the computation of all quantities reported in the examples discussed 
below. We stress that no approximations have been introduced into these equations. The 
simplifications we obtain are due to the fact that we only consider the infinite time limit 
in eq. 40. The numerical advantages of eq. 40 relative to direct matrix exponentiation 
are many fold. First, it is not necessary to pick a time to evaluate your expressions and 
somehow confirm that this time is both large enough to insure the steady state yet small 
enough to avoid numerical instabilities. Eq. 40 assumes t — > oo. Using this method one 
only has to find the eigenvalues and eigenvectors of the matrix L for a given excitation 
frequency to obtain both the intensity and the Q parameter. This matrix is three times 
smaller in linear dimension than the matrix that must be exponentiated to solve eq. 31. 
If higher moments are required, you still have only to diagonalize the L matrix for use in 
expressions similar to eq. 39. Finally, while matrix exponentiation requires that you repeat 
the entire calculation to obtain statistics for various detection possibilities (broadband, a 
single transition counted, several transitions counted, etc.), the present scheme only requires 
a single diagonalization for all possible detection schemes. Different detection possibilities 




(40) 
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2 Ep.E.U0^ +T ^ +r ps, 
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manifest themselves only through the matrix L + which does not have to be diagonalized. 
The pieces of eq. 40 dependent on matrix diagonalization (X ,Uo,p s . Sm ) do not vary with 
different detection schemes. This is a significant computational advantage when calculating 
emission spectra since the bulk of the calculation need only be performed a single time. 

III. CHROMOPHORE COUPLED TO A TWO LEVEL SYSTEM 
A. Model description 

As a first example, we consider the case of a chromophore coupled to a two level system 
(TLS). The two level system model is of interest both for theoretical reasons (it is arguably 
the simplest case of dynamics beyond that of an isolated two level chromophore) and also for 
its utility in describing the thermal behavior of low temperature glasses [52, 53]. The model 
is also frequently applied to the spectroscopy of chromophores embedded in low temperature 
glasses [54] . Although TLS dynamics is often treated as a purely stochastic perturbation of 
the chromophore system, we adopt a more precise, quantum mechanical picture here. The 
following description of coupled chromophore-TLS dynamics is quite terse. We refer readers 
to the review by Silbey [54] for more detail on the Redfield dynamics that we employ. 

The nature of TLS dynamics within the glass is presumably the localized rearrangement 
of a small cluster of atoms [52, 53] corresponding to movement between two distinct energy 
minima. The coupling between TLS and chromophore enters as a different effective split- 
ting between chromophore ground and excited states depending upon which minima the 
TLS resides in. Assuming this coupling is due to strain dipole interactions between chro- 
mophore and TLS we expect the interaction to scale as 1/r 3 in the distance between TLS 
and chromophore centers [54] . The basis of TLS "minima" states is not expected to be di- 
agonal as tunneling may occur between minima. In addition, coupling between the TLS and 
long wavelength phonons in the glass acts as mechanism for coupling the TLS-chromophore 
system to its glassy environment. Adopting the notation of sec. II C The mathematical 
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TLS 

FIG. 1: Energy level diagram for the composite chromophore -TLS system 
formulation of this picture is [54, 55] 

HuJ eg (A OL \ TLS , J TLS { A~\\ 

H « = -— + \2-4pi) a * + 2 a * ' (41) 

H e = + ^ + (^ + ^y™ + j_ a TLs 

q 

H b = 52blb q hu q . 

i 

Here, A and J are respectively the asymmetry and tunneling matrix element for the TLS 
and crJ LS ,&x LS are Pauh matrices in the basis of TLS localized "minima" states. u eg is 
the chromophore transition frequency in the absence of interactions. The index q labels the 
phonon modes of the system and b q , b q , u q and g q are the the creation operator, annihilation 
operator, frequency and TLS strain field coupling constants for the qth mode. 

We diagonalize the chromophore-TLS portion of our Hamiltonian and label the four 
eigenstates \a), \b) \c) and \d) (see fig. 1) in order of increasing energy (we assume fkj eg to 
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be by far the largest energy scale in the problem). In this basis eq.41 can be written 



H 9 = 
H e = 
V = 



uj a \a)(a\ + u b \b)(b\, 
w c \c){c\ + LO d \d){d\ 



( 6 -*+ 6 *) 

q 

-(\c)(d\ + \d)(c\) 



J 



(\b)(a\ + \a)(b\) + 



where cu a , u>b, u> c , u>d, oj 9 and u e are the frequencies 

= -\^ g + \^J 2 + {A-P)\ 
uj c = + 1 -u; eg -^Ji + (A + Py, 
Wd = +^ e , + iv/J 2 + (A + P) 2 

— — 



(42) 



(43) 



where we have set y = ^3. Note that we have intentionally omitted all (system) diagonal 
contributions to the system-bath coupling since these terms will yield no contribution to the 
Redfield matrix. 

Specification of H. is quite simple (if tedious) and proceeds by calculating the terms 
specified in eqs. 25 and 26. Since the bath is formed by a set of bosons (phonons), evaluation 
of the correlation functions is dictated by the well known properties of these operators. In 
particular since 



b q (t) 
bl(t) 



£ -iU)gt 



b q (0) 
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the correlation functions become 



(45) 



The coupling constants g q are chosen to reflect strain field coupling between TLS and the 
phonon bath [54]; they scale with q as q 1 ^ 2 . The ij and kl suffixes on g q indicate that there 
are additional constants that need to be included - either J/u e or J/u g depending upon 
which specific terms the indices refer to. Integration in time over these terms as specified by 
eq. 26 serves to create a delta function in frequency which makes evaluation of the sum over 
q trivially easy if we approximate the sum as an integral. By this approach we calculate, for 
example, 



■i;bb — e^^bb; 
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(46) 
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CO, 



L0 r , 



Where C is a collection of constants incorporating the coupling strength between TLS and 
bath, which is typically taken as a parameter used to fit experiment rather than estimated 
from first principles [14]. Of course the top two lines just express the phonon assisted 
transition rates from state d to c and b to a as expected. Other elements follow similarly. 
We make no effort to implement the customary secular approximations to these equations as 
the equations are solved numerically and the highly oscillatory terms will remove themselves 
from consideration naturally. 



B. Numerical results 



In this section we present numerical results for the model system described above. The 
framework for calculating the fully quantum dynamical results are spelled out in sec. II. 
Physical constants have been chosen to correspond with typical situations for a glassy ma- 
terial [14, 55]. In order to compare with our previous work on stochastic models, it is 
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necessary to map the above quantum description to a stochastic picture. Details for calcu- 
lating photon statistics for a stochastic TLS coupled to a chromophore has been presented 
in detail elsewhere [37]. Readers are referred there for a discussion, where we have employed 
notation identical to the present work. Determination of appropriate model parameters for 
the stochastic model, based upon the above quantum picture, is well established [54]. In 
the stochastic picture the TLS acts solely to modulate the transition frequency of the chro- 
mophore, causing hops between u eg + v and uj eg — v. The rate of hopping is given by for 
transitions to the less thermally occupied TLS state and R± for the reverse direction. The 
difference in energy of the two TLS states is provided by detailed balance. Correspondence 
with the quantum model is accomplished by 



The idea of the stochastic approach is that coupling between TLS and chromophore 
only manifests itself through modulation of the absorption frequency of the chromophore as 
modulated by TLS hops. TLS dynamics and thermal properties are completely unaffected 
by the chromophore, hence the total independence of TLS energy scale and flip rates on 
chromophore properties - i.e. these quantities are calculated by setting the TLS-chromophore 
coupling constant a to zero in our earlier expressions. Of course it is crucial to keep a in 
the frequencies, otherwise the TLS would have no effect on the chromophore at all. The 
stochastic approximation is expected to work quite well when a is small, in that case 
transition elements of the Redfield matrix are well approximated by using rates inferred 
from eq. 47. It should be noted that the stochastic approach is obviously deficient in one 
sense. There are four possible transition frequencies implied by the quantum level diagram 
in fig. 1 and the stochastic picture only predicts two. For small a and/or large r, half the 
transitions rarely occur because of poor Franck-Condon overlap. Given our notation, the 
transitions c — > a and d b are the strong ones (assuming weak coupling). At high couplings 
strengths, half of the transitions will necessarily be missed by the stochastic picture. The 
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following numerical examples highlight both the practicality of the present fully quantum 
approach in calculations as well as the shortcomings of the popular stochastic approximation 
over certain parameter regimes. 

1. Weak coupling between chromophore and TLS. 

"Weak" coupling between the chromophore and TLS is dictated by the condition A ^> 
^§ = P. Physically this can result from either a small coupling constant a or a large 
distance between the chromophore and TLS. As discussed above, in this case, results of the 
quantum model and stochastic model should be quite similar (at least for the line shapes 
[54]). In the left panes of fig. 2, we present the long time lineshape and corresponding 
Q parameter spectrum for the case of slow TLS modulation and weak TLS-chromophore 
coupling. The physical constants chosen are detailed in the figure caption and represent 
realistic numbers for an organic dye molecule embedded in an amorphous host [55]. We 
compare the quantum model with the associated stochastic approximation. As expected, 
the line shapes for the two approaches are identical at the resolution of the figure. The two 
peaks represent the two optical transitions with appreciable overlap (a — > c and b — > d). 
The other transitions are so weak as to be invisible at this scale. The difference in peak 
heights is due to the difference in thermal occupation probabilities for the two TLS states 
(which are basically unmodified by chromophore state due to the small value of P in the 
quantum model). Peak shape is Lorentzian with both linewidths given by the spontaneous 
emission rate (full width at half maximum is r ). The TLS flipping is so slow in this case 
that it contributes negligibly to the linewidths. 

The right panes of fig. 2 display similar information to the left, but with parameters 
chosen to insure that the TLS flip rate is faster than the difference in transition frequencies, 
v. For simplicity we increased the flip rate by increasing the value of C. While this is 
physically questionable, it does provide the only direct means to increase the TLS flip rate 
while leaving all other behavior identical. In this case, the lineshape consists of only a single 
peak due to motional narrowing of the optical transition [54, 56]. As in the slow modulation 
limit, we find quantitative correspondence between stochastic and quantum models for the 
lineshape calculation. The stochastic model does deviate slightly from the quantum result 
in the calculation of the Q parameter. Though the deviation is slight, it is interesting to 
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FIG. 2: Absorption line shapes and Mandel's Q parameter spectrum in the limit of weak coupling 
between the chromophore and TLS. Lineshapes are presented in arbitrary units. Left and right 
halves correspond to slow and fast modulations respectively. Physical parameters used in this 
calculation include Tq = lOOMs -1 , f^o = lMs~ , T = 1.7K and quantum model parameters taken 
from ref.[55], namely, A = 2.8K, a = 3.75 x 10 11 nm 3 s _1 , r = 5.72nm, J = 3 x 1Q~ A K. For the slow 
modulation we used C = 3.9 x 10 8 s~ 1 K" 3 while for the fast modulation C = 3.9 x lO^s" 1 ^ -3 . 
Within the stochastic approximation these numbers translate to (eq. 47) v = 1.02 x 10 9 s _1 and 
E = 2.8K. The upward flip rate R-\ = 23.5s -1 for the slow modulation and 2.35 x 10 11 s _1 for 
the fast modulation. In the slow modulation, no discrepancy between quantum and stochastic 
treatments is found. For the fast modulation the line shape is the same for both quantum and 
stochastic treatments while in the Q parameter there is a small difference between the models. The 
inset focuses on this difference. 



note that there are cases where the stochastic model is perfect for lineshapes, yet imperfect 
for higher order statistics. All in all though, for weak coupling, the stochastic approximation 
is seen to perform well both at slow and fast TLS modulation rates. We note that in the 
limiting cases of slow and fast modulation displayed here, the observed spectra can also be 
predicted on the basis of the physical approximations introduced in ref. [39]. 
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2. Strong coupling between chromophore and TLS. 

"Strong" coupling is insured by the condition A ~ P — ^S' - this case, the quantum 
model differs from its associated stochastic approximation in both line shape and Mandel's 
Q parameter. The left panes of fig. 3 display results for the strong coupling and slow 
modulation parameter regime of both the quantum and stochastic dynamic treatments. 
In contrast to our earlier example, strong coupling now implies that transitions between 
states d — > a and c — > b are important and occur with some finite probability within the 
fully quantum treatment. Since peak widths are smaller than interpeak spacing, peaks 
corresponding to all four possible transitions are clearly visible in the quantum mechanical 
modeling. The relative height of the two central peaks in the line shape are (as in the 
previous example) related to TLS thermal occupation probabilities. Since E <C kT for 
the chosen parameters, both central peaks have effectively the same height. The intensity 
of the outer two peaks is predicted based on the probability to excite an "off diagonal" 
transition (a — > d or b — > c) relative to diagonal transitions. Mathematically this probability 
is dictated by the square of the Rabi frequency for the transition in question. Equivalently 
(see eqs. 21 and 22), the ratio of the left two peaks or the right two peaks is predicted to 
be T db /T da (1.94 for the case shown), which agrees with the numerical results. It is obvious 
that the stochastic approximation predicts a very different line shape and Q parameter since 
it doesn't account for the transitions d — > a and c — > b. While one could argue that the 
stochastic model does do a good job in predicting that portion of the absorption lineshape 
which it is capable of reproducing (the center two peaks), even the center two peaks are 
clearly off in magnitude for the Q parameter. The stochastic model fares very poorly in this 
parameter regime (strong coupling, slow modulation). 

The failure of the stochastic model in this case was predictable and we can trace its origins 
back to failures to reproduce the full system dynamics in a realistic manner. The right panes 
of fig. 3 are meant to display that we understand exactly where these failures are coming 
from. These panes actually display two different cases (although they overlap so only a 
single line is visible): the stochastic calculation from the left panes and a modified quantum 
calculation where the evolution operator was altered such that all non-diagonal transitions 
were turned off (Q ad = Q bc = T ad = T bc = and Q ac = Q bd = Q and T ac = T bd = T ) 
and all Redfield elements were calculated assuming that u g = u e = E. While these two 
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FIG. 3: Left panes: The line shape and Mandel's Q parameter spectrum for slow TLS modulation 
with strong coupling between the chromophore and TLS. Due to the strong coupling, d — * a and 
c — * b transitions are significant within a fully quantum framework and result in two additional 
peaks relative to weak coupling results. The stochastic approach completely misses these additional 
spectral lines and fares poorly in reproducing the magnitude of peaks in the Q spectrum. The plots 
correspond to quantum model parameters: Tq = 40Ms -1 , Qq = O.lMs , T = 1.7K, A = 0. 006.fr, 



J = 0.008^, C = 3.9 x 10 8 K~ 3 s-\ 



a 



3.75 x 10 n nmV 



5.72nm. Corresponding stochastic 



parameters are: v = 501Ms -1 , E = 0.01K and i? T = 42307s" 1 . The right panes display that it 
is possible to reduce the fully quantum mechanical treatment to the stochastic results by turning 
off half of the allowed transitions and calculating Redfield elements in a manner consistent with 
the stochastic approach (see text). In other words, it is relatively simple to trace the failures of 
stochastic modeling. 



changes do not fully reduce the quantum calculation to the stochastic treatment from a 
mathematical standpoint, the physical basis is clear. The alterations explicitly remove the 
non-diagonal transitions that the stochastic model necessarily misses and it evaluates the 
TLS jump rates in the same approximation inherent to the stochastic approach. There are 
more subtle effects within the Redfield treatment (as in the evolution of coherences) so that 
our ad hoc alterations do not fully limit to a stochastic model, however these effects clearly 
do not contribute to the lineshape and Q spectrum calculations. The primary problem with 
a stochastic model in predicting photon counting observables is in the loss of "off-diagonal" 
nuclear transitions and incorrect estimation of relaxation rates. 

In fig. 4 we show two cases of reasonably fast modulation speed and strong coupling; 
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FIG. 4: Line shape and Mandel's Q parameter for intermediate TLS modulation rate, with "strong" 
coupling between the chromophore and TLS. The quantum model parameters are the same as in 
fig. 3 except for the coupling constant which is modified to C = 3.9 x 10 12 K~ 5 s~ 1 , corresponding 
to upward flip rate i?| = 4.23 x 10 8 s~ 1 in the stochastic model. In the left panes the Rabi frequency 
coefficient is Qq = 10 5 s _1 , while in the right panes f2o = Tq = AOMs" 1 . Comparison of the left and 
right panes shows that antibunching increases as excitation and emission rates become comparable. 

the difference between left and right panes is quantitative (see the figure axes for Q) and is 
intended to display the fact that you can tune the Q parameter by adjusting field strengths. 
For a simple two level chromophore, antibunching is maximized when excitation and emission 
rates are equalized [57] and a similar effect is seen here. Although both quantum and 
stochastic models will eventually narrow into a single peak for high enough flip rates, it is 
interesting to see in this intermediate regime that the stochastic model has already narrowed, 
while the quantum picture retains a more complex structure. This structure is visible in 
both the lineshape and Q parameter calculations. 



3. Emission spectra 

In fig. 5 we display emission line shapes and Mandel's Q parameter spectra for the 
same physical parameters selected in fig. 3 (excepting the Rabi frequency, which was set 
to provide relatively large magnitudes of the Q parameter in the anti-bunching regimes). 
As discussed previously, our simulation methodology does not allow for true calculation of 
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emission spectra. The frequency dependence we obtain is resolved solely on the basis of 
individual state to state transitions - we assign all photons emitted for a given transition the 
resonance frequency of that transition. Hence, the "lineshapes" in fig. 5 are not broadened 
by the radiative lifetime of the chromophore or by any other source and line shifts are not 
captured. Physically, the spectra we obtain would match an experimental measurement with 
an instrument unable to resolve frequency differences less than the radiative line width. 

The multiple panels in both rows of fig. 5 reflect different laser exciting frequencies. Four 
different resonant excitations corresponding to all possible transitions and two off resonant 
frequencies are considered. Clearly, there is a strong dependence in the emission spectra 
on the exciting frequency. This is expected since TLS dynamics are slow enough in this 
problem that the TLS does typically not have a chance to relax to equilibrium while the 
chromophore is excited. Resonant excitation to state c, regardless of which ground state 
(a or b) the transition starts from results in the same emission line shape (left two panes 
of the top row of fig. 5). The relative peak heights simply reflect Condon overlaps in the 
spontaneous emission process from state c back to a or b. These overlaps don't care how state 
c was excited and generate identical emission spectra regardless of which resonant transition 
is excited. Similar arguments explain the right three panes of the top row of fig. 5. All 
three excitation frequencies result in the occupation of state d and the emission lineshapes 
are insensitive to details of the excitation beyond this fact - even when the excitation is 
off resonance with either a — > d or b — > d transitions. When an off resonant excitation 
is considered that has equal probability to excite to either c or d, the emission lineshapes 
reflect a symmetric combination of the previously discussed cases (third pane of the top row 
of the figure). 

In contrast to the lineshapes, Q parameter spectra are highly sensitive to excitation 
frequency (bottom row of fig. 5). The basis for this effect is quite simple. When photons 
are counted at the same frequency of the exciting laser we expect to see photon bunching. 
For example, looking at the leftmost peak in the leftmost pane of the bottom row we excite 
b — > c transitions and monitor c — > b emissions. Photons are repeatedly ejected as this 
cycle repeats until spontaneous emission induces a c — > a transition (or the TLS flips), at 
which point the system is off resonance and has to wait for a TLS flip to return the system 
to the excitable state b. The interspersion of bright and dark intervals leads to bunching 
phenomena and a positive Q parameter. In contrast, when excitation does not correspond 
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FIG. 5: Emission lineshapes and the Mandel's parameter Q for slow modulation limit with "strong" 
coupling between the chromophore and TLS. The excitation laser frequencies are marked in the 
figure using "j" . The excitation frequencies, from left to right, are u> eg +uJcb, u eg +oJ C a, u eg , oJ e g+^db-, 
uj eg + LJda, and uj eg + 0.6u>d a (see fig. 1). The spontaneous emission rate and the Rabi frequency 
are Tq = 40Ms _1 , Qq = AMs , respectively. The quantum model parameters are: T = 1.7 K, 



A = 0.006K, J = 0.008, q = 3.75 x 10 n nm 



"\ C = 3.9 x lO 8 ^ 3 ,^ 1 and r 



5.72nm. 



to the monitored transition (second peak from left in the leftmost pane) a three state cycle 
repeatedly occurs (b ^ c — > a or a similar variant) as photons are detected. There 

is no jumping between periods of "bright" or "dark" since the pathway for repeated photon 
emission necessarily involves both TLS and radiative/excitation dynamics. The chosen 
timescales in this example insure that no single rate is limiting over all others in this cycling 
process and antibunching results (if a single timescale were completely dominant we would 
expect Q — 0). Similar arguments can be applied to the remaining panes of the Q parameter 
spectrum. This example makes a clear case for measurement of higher order photon counting 
moments. Different aspects of system dynamics are captured in the measurement of the Q 
parameter beyond what is seen in simple lineshape statistics. Furthermore, examination of 
the emission statistics provides a more detailed measure than possible solely on the basis of 
absorption statistics. 
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IV. A CHROMOPHORE WITH NUCLEAR VIBRATIONS COUPLED TO AN 
HARMONIC BATH 



A. Model description 

As a more complex example of multilevel quantum dynamics we consider the case of a 
chromophore with an harmonic vibrational degree of freedom. Coupled to this vibrational 
coordinate is a bath modeled by an ensemble of harmonic oscillators. Such models are 
standard in the treatment of molecular spectroscopy [47], but have seen little prior use 
in the treatment of photon statistics. Within the Born-Oppenheimer approximation, the 
Hamiltonians of the chromophore in its electronic ground, \g), and excited, |e) , states are 
taken to be 

H g = ^hw [P 2 + X% (48) 
H e = hw eg + ^fko [P 2 + (X -X ) 2 ], 

where X and P are related to the nuclear position coordinate x and momentum p by 



mu h 

The vibrational coordinate thus has frequency u and ftw eg is the excitation energy for the 
— transition. x is the shift in equilibrium position of the nuclear coordinate between 
excited and ground states (see fig. 6). The interaction with the thermal bath is assumed to 
be linear in both X and bath coordinates Xi, i.e., 



V = CXY,*i (50) 



3 



where C is a constant specifying the interaction strength between system and bath. The 
harmonic bath Hamiltonian is 

H B = Y l %i{ P i+ x j)- ( 51 ) 
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The above definitions of H g , H e , V and Hb provide all necessary information to proceed 
directly with the calculation of L and related quantities as detailed in section II. We make a 
few brief comments related to the calculation of Redfield elements below in order to clarify 
our notation. More detailed presentations can be found elsewhere [47, 48, 58]. 

The linear interaction between bath and system in only capable of effecting transitions 
between adjacent vibrational states in the same electronic manifold, i.e. \n) — > \n + 1) or 
\n) — > \n — 1). This is seen, by introducing the usual creation and annihilation operators 
(a = (X + iP) I a/2, a* = (X — iP) / \/2) to write the interaction matrix elements between 
excited state levels in the form 

Vn e ,n' e = ^C[Vrf e Sn e ,n> e -l + \/^i,»' e +ll J^( a j + ( 52 ) 

3 

and similarly for the ground state. The creation and annihilation operators only allow for 
adjacent transitions as indicated by the above delta functions. The bath properties 

aj (t) = e-^S(O) (53) 
a ){t) = e^'aJ(O) 

( aj a]) b = {l-e-^y 1 
(a] aj ) b = e-^il-e-^y 1 . 



are used to evaluate all correlation functions associated with the Redfield matrix calculation. 
In this model the interaction matrix V is explicitly real leading to a slightly simplified 
calculation for the Redfield matrix 

^mnpq tpmnq tqnmp $rnp ^ ^ ^nrrq $nq ^ ^ t mrrpi (^4) 

r r 

where 

1 f°° 

tpmnq = ^ J (V pm (r)V nq (0) ) b . (55) 

tpmnq is non-zero only if both of the pairs (p, m) and (n, q) involve states in the same electronic 
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FIG. 6: Schematic description of a system with two electronic levels and single harmonic vibrational 
mode. 



manifold. The integration can be carried out and yields 

_ r i ~ r sr^ 5(ujj + u qn ) + e- /3huJ ^5(uj j - u; qn ) 

tpmnq ~ -^0 [V ?™>p,m-l + \fPO P) m,+l\ [V nO q>n -i + qO q , n+ i\ ^ 1 _ -fffaj ■ ' 

(56) 

Unlike the TLS model, in this case every allowable u qn is exactly the same and is equal to 
c^o- This is due to the equality of spacing between levels in the harmonic oscillator model 
and the form of V which only allows for adjacent transitions. Thus, the density of bath 
states is not important in calculating the Redfield matrix elements in this case and only 
a single constant Rq enters into the Redfield description as a measure of coupling between 
system and bath. For example, elements of the form R nn;n+ i n+ i are given in our notation 
by 

Knn;n+ln+l = 2i? (™ + 1)(1 " e"^)" 1 . (57) 

Since this element reflects the rate of transition from harmonic oscillator state \n+ 1) to \n), 
it is clear that i? is closely related to the relaxation rate of our vibrational coordinate. 



B. Numerical results 



In the following calculations we choose physical parameters specifying the chromophore to 
be c^o = 3.77x 10 13 s _1 , Xq = O.llA, m = 10 5 m e (m e is the electron mass) and T = 10K (The 
energy difference between neighboring levels of the harmonic oscillator Hujq corresponds to 
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FIG. 7: The line shape and the Mandel's Q parameter spectrum as a function of exciting laser 
frequency for a chromophore with an harmonic vibrational coordinate. The spontaneous emission 
rate and Rabi frequency are To = = 10 8 s _1 and the coupling strength is Ro = 10 7 . Physical 
parameters specific to the chromophore are detailed in the text. 

temperature of 287. 75K). While these numbers are suggestive of a heavy diatomic molecule 
(like I2) in a low temperature matrix we have not made a serious attempt to connect these 
calculations with physical systems. Rather, we have chosen x to provide Condon overlaps 
that are close to vertical, while still insuring finite probability for transitions up to — 6. 
We have also set the temperature somewhat arbitrarily while we will freely adjust Rq in 
the following examples to meet our needs in displaying various phenomena. The Redfield 
approach we employ is necessarily limited to a finite number of states due to numerical 
considerations. We can not solve the equations for N = 00. In the numerics presented here 
we used 10 levels in each of the electronic states (n = to n = 9). It was verified that 
altering the number of vibrational states to include more levels did not change any results 
at the resolution of the presented figures. We note that the size of L for these calculations 
is 400 x 400. Using the methods of sec. II C requires only diagonalization of this matrix, 
which is a simple task for modern computers. 
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1. Weak coupling between system and bath (Ro small) case. 

The case of weak coupling corresponds to slow vibrational relaxation. In fig. 7 we show 
the line shape and the Q parameter for a case in which the relaxation rate is slower than 
all other rates in the problem including the spontaneous emission rate, Rabi frequency and 
oscillator frequency. This leads to non-thermal distributions of vibrational levels within 
both electronic manifolds at steady state since the system is unable to fully relax between 
subsequent photon emission/absorption events. Interestingly, the variation of these steady 
states with excitation frequency and the variation of Condon overlaps between the various 
transitions leads to Q parameter values spanning a range of positive and negative values 
depending upon the excitation frequency. It should be noted that although the spectra 
appear to have only been evaluated at the various allowed resonance frequencies, this is not 
the case. It is simply the case that the radiative linewidths are very much narrower than 
discernible at the resolution of the figure. 

Fig. 8 shows the line shape and Q parameter for a case in which the relaxation is slow 
relative to the harmonic oscillation frequency u , but is faster than the spontaneous emission 
rate and the Rabi frequency. In this case the relative amount of power absorbed by each 
possible transition is expected to agree with linear response predictions since the vibrational 
state of the chromophore should almost always be in the relaxed (n = 0) state without 
significant perturbation by the relatively weak coupling to the field. Linear response theory 
predicts that the strength of each transition is due to the Condon overlap between n = in 
the ground state (remember kT <C hu in this model) and the various excited states. The 
displayed lineshapes appear to contradict this prediction, most clearly due to the very tall 
zero phonon peak at ojl = oj eg relative to the other peaks. However, the height of this line 
is due to the fact that this transition is not broadened by non-radiative processes as are the 
remaining transitions. The linewidth of the — line is approximately equal to T whereas 
the other widths are dominated by non-radiative decay on the order of Ro and are 100 times 
wider. The relevant quantities to compare with linear response results are the intensities of 
each transition integrated over the local vicinity of the transition. In fig. 10 we display such 
integrated absorption peaks alongside emission lines (discussed below). These integrated 
lines show perfect agreement with linear response results with relative intensities directly 
proportional to the square of nuclear overlap. 
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FIG. 8: Similar to fig. 7, but with R = 10 10 s _1 , T = W 8 s~ 1 and tt = 10 6 . This system is in the 
linear response regime. The inset shows the variation of Q in the vicinity of ujl = w eg . 

One interesting point to note about the Q parameter in these calculations is that it 
undergoes rapid variation with excitation frequency in the vicinity of the — line. While 
this behavior does not seem amenable to simple explanation, it has been observed previously 
in simpler models both numerically [38] and analytically [42]. It should also be emphasized 
that the magnitude of Q is largely due to the ratio between T and Q as seen in fig. 4. 
Here this ratio is large, leading to small negative Q values. Smaller ratios lead to larger 
magnitudes of Q (when Q is negative). 



2. Strong coupling between the system bath (Rq not small) case: 

An example of fast relaxation, with R on the order of Uq, is shown in fig. 9. In this 
case the width of the peaks is of the same order as the distance between the peaks and line 
shape is clearly not a series of thin sticks as in previous examples. Note that since the peak 
at lol = uo eg does not involve any thermal relaxation it is independent of Rq. The width of 
this peak is still specified by the spontaneous emission rate, which is orders of magnitude 
lower than the remaining peak widths (on the order of Rq) leading to its very large height. 
In this plot we have chosen identical values for Tq and Qq, which leads to sizable negative 
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FIG. 9: Similar to fig. 7, but with To = = 10 9 s _1 and Rq = O.Iloq (solid line) and Ro = 0.05<jJo 
(dashed line). Note that in this case the width of the peak at is much smaller than the width 
of the other peaks, since it does not depend on Rq. Peak widths are given by the non-radiative 
lifetime of the various states for all other transitions. 



Q values for the — line. 



3. Emission spectroscopy 

In fig. 10 we show the emission line shape and Q parameter spectra for parameters 
appropriate to the linear response regime (identical parameters to fig. 8). It is shown that 
in this case the line shape is the same for all excitation frequencies (in the fig we show 
uj l — u) eg = 0, 2ou , Aujq). It is also shown that integration of the absorption spectrum over 
the individual transition linewidths provides a mirror image of the emission line shape as 
expected in the linear response regime. Recall that our emission line shapes are sensitive only 
to individual transitions, so the emission spectra are automatically of the "integrated" type 
and comparison between emission and integrated absorption is completely natural. While 
emission lineshapes are insensitive to excitation frequency in this regime, the Q parameter 
exhibits strong dependence on excitation frequency. 

Fig. 11 shows the emission line shape and Q parameter for stronger driving fields and 
slower relaxation rates than present in fig. 10. The system is no longer in the linear response 
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FIG. 10: The emission line shape and emission Q spectrum for parameters reflecting the linear 
response regime. Three different excitation frequencies are considered as noted in the legend. The 
chosen physical parameters parallel those of 8. Since the system behaves in accord with linear 
response, the emission line shapes are the same for all excitation frequencies and also in agreement 
(mirror image) with the integrated absorption spectrum. 

regime and line shapes differ for different excitation frequencies. The parameters were chosen 
to equalize all relevant physical timescales, demonstrating that there is no simple relationship 
possible between excitation frequency, emission lineshape and emission Q parameter possible 
in general. 

V. CONCLUSION 



We have introduced a practical framework for the calculation of photon counting statistics 
in quantum systems with multiple levels and dissipative coupling to a thermal environment. 
The present scheme generalizes previous work by extending the treatment of chromophore 
dynamics beyond the stochastic models historically applied to single molecule spectroscopy. 
Our model calculations for TLS dynamics explicitly demonstrate some of the failings of 
traditional stochastic modeling. In the case of harmonic vibrations, use of a stochastic 
model is even more suspect since all quantization of the vibrational coordinate will be lost. 
Although one could envision more elaborate kinetic schemes in an attempt to model these 
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FIG. 11: Emission line shape and Q spectrum for parameters outside the linear response limit. 
i?o = To = S7o = 10 8 s _1 . The exciting fields are as in fig. 10. In this case the linear response 
approximation is not valid anymore and both the line shape and Q spectra differ with excitation 
frequency. 



systems, it seems more straightforward to simply treat the dynamics correctly, quantum 
mechanically, from the outset. The methods presented here provide a prescription to do 
this. 

We acknowledge that there is an unfortunate amount of machinery behind the calculations 
that we have presented here, however it is important to stress that 90% of this machinery 
is associated with the implementation of the Redfield formalism (calculation of the matrix 
L in our notation). Eqs. 40 are very simply applied once L is given; simply diagonalize 
the matrix and perform a few simple matrix multiplications as implied by the formulae. 
The generating function approach, while necessarily encumbered by the usual difficulties in 
simulating dissipative quantum systems, adds no new significant conceptual or numerical 
problems. Photon counting statistics are therefore readily available at no more expense than 
normally expected for calculation of density matrix dynamics. This remarkable fact seems 
to be the strongest point in support of the generating function methodology. 

Several of our calculations have presented results for emission spectra and the corre- 
sponding Q parameter quantities. Although such measurements are not yet within the 
capabilities of experiment, we believe that a strong case can be made for the development 
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of single molecule detectors with spectral resolution. It is clear from our model calcula- 
tions that emission spectroscopy provides a different and (when combined with absorption 
spectroscopy) more revealing signature of chromophore dynamics than obtainable from ab- 
sorption alone. This is not surprising, but the present study is (to our knowledge) the first 
to demonstrate this fact explicitly. As we have repeatedly stated, the present scheme for 
emission spectroscopy is sensitive only to molecular transitions and not directly to emis- 
sion frequency. Emission frequency is assumed to be on resonance with specific transitions. 
While this approach works well in the limit of weak coupling to the environment, stronger 
coupling invariably leads to level shifts, motional narrowing as associated complications. A 
general and practical formulation of true emission photon counting statistics has yet to be 
developed. 
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